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ABSTRACT 


To aid in the modeling of a steady state spin, the 
equations of motion of an airplane are formulated ina 
cylindrical coordinate reference frame. The derivation 
Sr the equations 1s presented and the resulting equa- 
tions are simplified for the equilibrium spin condition. 
These Simplified equations are used in an unconstrained 
computer parameter optimization technigue that algebrai- 
cally solves the differential equations for the equilib- 
rium state. The results of the computer work are 
presented and compared with previous prediction schemes. 
The potential of the method is demonstrated by applica- 


tion to a study of the effects of density variation. 
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Prot eOhaoYMBOLS 


The definitions of the symbols used in this paper are as 


follows: 


WiisiigiesS Dall, 96n c 
Mean aerodynamic chord, ft 


Lagrange generalized force, lbs or 
ft-lbs 


Aerodynamic forces along the body axes, 
lbs 


Pe@eleracivon Of Gravley, EW eee 


Body axes moments 
inertia, slugs-ft 


,and PacCauctus OF 


Rolling moment acting about the X 
body axis, ft-lbs 


Pitching moment acting about the Y 
body axis, ft-lbs 


Yawing moment acting about the Z 
body axis, ft-lbs 


Airplane mass, slugs 


Angular rates of rotation about the 
X, Y and Z body axes, rad/sec 


Dynamic pressure, 1/2 a 
Lagrange generalized coordinate 
RaCiuUSswOpenerrcal path of aircraft, £t 
Wing area, Bea 
Lagrange kinetic energy, ft-lbs 
imeruial velocity of the C.G., ft/sec 
Eageange eoOtential energy, ft—lbs 


Aircraft body axes 





Ky, 1S) 


p 
W,O49 
(") 


CarteSian coordinate reference axes 
that are attached to the cylindrical 
maGdmMus vector. 

MWeerarie altitude, ft 

Angle of attack, degrees 

Direction cosine, dimensionless 
Angle of Sideslip, degrees 


Cylindrical orientation angle, radians 


Aileron deflection (positive when right 
aileron trailing edge is down), degrees 


Elevator deflection (positive when 
elevator trailing edge is down), degrees 


Rudder deflection (positive when rudder 
trailing edge is left when viewed from 
above), degrees 

Atmospheric density, Sinceyee 


Euler orientation angles 


Derivative with respect to time 


Geeteterents and Gerivatives: 








Cc, = "x 
ma qo) 
Cc ne 
de 36. 
a 
Gc eg 
ds 90. 
ae 
P _ dC y 
X pb 
p 3 (5u7) 
2 - 9C > 
oe ea 
2V 

















a ee me 7 ng 
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e a 
0C 4 oC. 
C == Cc =) eee 
m qe n 5 90 
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aC 
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ieee ENT RODUCTION 


A. BACKGROUND 

As a result of post-—stall/spin accidents a Significant 
number of lost fighter aircraft and pilot fatalities have 
been recorded over the last decade. Investigations have 
shown that these aircraft exhibit poor spin characteristics 
and that recovery from a fully-—developed spin is usually 
difficult or impossible. For these reasons the spin is no 
longer considered a Paemeeat. iancuver and is regarded by 
pilots as an undesirable and potentially dangerous flight 
Gondition. 

Aemougm: the £Llight condations leading tO a Spin can 
generally be avoided, the combat pilot often finds the need 
to utilize the edge of the flight envelope in a combat 
Situation. More often than not, it is only in._this boundary 
area nate the combat superiority of one aircraft over an— 
other is realized. If the aircraft demonstrates unsatis— 
factory stall and spin characteristics, the pilot is 
POeeeoiined) or restricted from utilizing this outermost 
region of flight to his best advantage. As a result, the 
Overall weapon system effectiveness and margin of superi- 
Simibey son themalrcrart 1s lost. 

In light of the poor stall and Spin characteristics 
of current fighter configurations a more detailed consid- 


eration must be given to these characteristics during the 
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Saehyeadesign stages of future fighter aircraft. In recog- 
nition of this fact this research was undertaken to provide 
eem@re accurate analytic tool with which the designer can 
investigate the spin characteristics of future aircraft and 
to aid in understanding how various factors in design may 
affect these spin motions. It was not proposed that this 
effort solve the spin problem itself, but it was hoped that 
it would prove to be a useful approach and that it would 
provide information complementary to that obtained experi- 


mentally. 


Peo COPE 

Traditionally, analytic spin prediction has involved 
Parmer an approximated solution (by requiring that only 
part of the equilibrium conditions be satisfied) or by com— 
Pmeer generation of time-histories of characteristic para-— 
meters (by integrating the equations of motion forward in 
time nae the average values of the variables are approxi- 
mieelyecOnstant). More recently the idea of mathematically 
casting the equations of motion into a form specifically 
intended for spin prediction has shown some promise. Of 
particular note is the work carried on by Buehler (Ref. l) 
and Champoux (Ref. 2). These earlier works provided much 
of the initial groundwork for this research. 

In developing the analytic approach to the spin problem 
a set of objectives was established to subdivide the re- 


search problem and also to serve as guidelines to ensure 


11 


mietie tne work proceeded in the proper direction as effec— 


tively as possible. These objectives were: 

ies TOorsuitably model the aircraft in a spin. 

Peelomcdevelope the Spin equations of motion. 

3. To establish a suitable technique of solving the 
equations of motion for the equilibrium spin con— 
Gan teee@ ia, 

4. To develop a computer program that will predict 
areca temsoln CharacteriStics. In addition, the 
BeogmanmesniOould pe flexible enough to aid in the 
understanding of how various parameters effect the 
developed spin motions. 

5S. EO demonstrate the usefulness of the method and to 


verify the results with other prediction schemes. 
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Zi Dove LOPMENT OF THE ALRCRAFT EQUATIONS OF MOTION 


Pa COORDINATE SYSTEM 

To an observer on earth, the aircraft's center of grav—- 
ity in a Steady spin appears to trace a descending constant 
radius helical path about an imaginary vertical central 
axis. By asSuming that the rotation of the earth is negli- 
gible, an inertial reference frame can be established that 
is fixed to the earth's surface at the point of intersection 
eens imaginary Spin axis. 

The most convenient method of describing this motion is 
moeoughn the use of cylindrical coordinates. In this formu- 
lation, the central axis of the cylinder 1S superimposed 
on the imaginary aircraft Spin axis as shown in Figure l. 
The advantage of this choice for an inertial coordinate sys-— 
monmis tMat the center of gravity of the aircraft 1s con- 
veniently located in terms of an altitude coordinate (Zo), 

a radial coordinate (R) and an angular position coordinate 
(y). 

The body axes, as shown in Figure 2, were chosen for 
the formulation of the equations of motion. This allows 
the direct use of standard wind tunnel data which was avail- 
Helomter tmissresearch effort. Symmetry of the aircraft 
PeotlemEnom "and 2eaxes was assumed So that all cross—product 
Spmineneaia. terms, wlten Che exception of Loo! were identi- 


Galuliy Zero: 


age 








Pugtiaem!.) inertial Cylindrical Coordinate 
system. 





Bagure 22@Aircraft Axes. 
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The Orientation of the aircraft is described by intro- 
ducing another coordinate system and a set of three rota-— 
tion angles. A fixed carteSian coordinate reference frame 
(X,Y, 724) iPSoePpesittlonea at the center of gravity of the 
aircraft and aligned with respect to the inertial cylindri- 
cal frame so that the Za axlS 1S parallel with the vertical 
central axis and the Yy axis coincides with the --R position 
vector. Although this cartesian coordinate reference frame 
remains fixed in orientation with respect to the inertial 
cylindrical frame, it is still free to move in space as Zoe 
Rand y take on different values. The orientation of the 
body axes relative to the cylindrical frame of reference is 
geeven by the aircraft Euler angles (,6,¢), which define an 
ordered series of three consecutive rotations as shown in 
Figure 3. The three position coordinates (ZorRr¥) together 
with the three orientation angles (W,8,¢), allow a Bones 
@eseraption Of the motion of the aircraft in space and will 
serve as the generalized coordinates in the Lagrangian 


development of the aircraft equations of motion. 


B. AXES TRANSFORMATION AND AIRCRAFT ANGULAR RATES 

Before proceeding with the development of the equations 
of motion, certain algebraic relations are needed that will 
transform the actual forces, moments, velocities and angular 
rates of the aircraft, about its body axes, to the cylindri-— 
cal reference frame. These Pee coat ons are nothing more 


than vector projections from one system of axes onto another. 


i> 








INERTIAL 


Bagiuiee 3, ) Spin Model: 


~~ 


eonencee Enat “the body axes! are inzelally coincident 
with the eee aa eae (X)7Y)72)) and then are rotated 
about 24 by an angle ~ to form an intermediate reference 
frame (X5,Y¥5,25),5 as shown in Figure 4. 
Since both of these orthogonal co- 
ordinate systems form a basis for 
three-space, any point (P) can be 


described by a vector (R) from the 





Galeimieco that point in terms of 
Figure 4. wW Rotation. 


1G 





either system. 


Mathematically, this can be expressed: 


nw 
° 


r= X.i. + Y 7 ea 


MOI ots ~2'9 ~ 235 * 45% (1) 
but, 
em ae yD ty 445° Jy tt 444° Sy 
eee) ty ID Jy tT 24d” Sy fe) 
ea, = Xk,- 2, + vik: 5, + 2, k- ky 


which iS written more compactly in matrix form as: 


“Aw Aw “aw “Aw “w “Aw 


Xo do*t, 45°37 2 Q*ky a 
2 = pty) (Dee ID) Yi; ©@) 
2 eo eee oe oa 


The three—by—-three matrix of unit vectors 1S called a trans-— 


PeommiatlOn Matrix and denoted, in this case, by [T }] and 


a/m 
read "the transformation from frame one to two." GaieeyirnG 
out the indicated dot products for the rotation, the trans— 


formation matrix has the value: 


cos Ww oie 0 
[21] = sin v cos y 0 (4) 
0 0 1 


The second rotation, 8, about the Y. ax1lsS carries sub- 


scripted coordinate system two into system three, as 


7 





shown in Figure 3. The resulting transformation matrix 


nS : 


cos 8 0 Sin 6 
[73/2 | = 0 - il 0 (5) 
“-Sin 86 0 cos 8 


The final rotation, 6, about the X4 axis carries sub- 
Scripted coordinate system three into the unsubscripted 


body axes, aS shown in Figure 3. This resulting trans— 


menmation Matrix is: 


ik 0 0 
"p73 | == 0 cos @¢ Sin @¢ (6) 
0 i—-Sin cos ¢ 


These individual transformation matrices can be multiplied 
together ons a Single three—by-three matrix which will 
be denoted by (By). Since each transformation matrix 
premultiplies the vector arrived at in the previous step, 
the total eeencromiaion from the cartesian coordinate 


frame to the aircraft body axes is given by: 


al - 273] : 3/2] , Pa i's 
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Sieeyine Out the Matrix multiplication, 


formation has the following form in which Che 


elements: 


CO cos0cosy 


le 


Sindsind 
cosy 
—cos¢dSiny 


ey) = = |%21> 


3 4= cos¢sindg 
cosw 
+singsiny 


Cy o> cos0sinw 


W595 Singsind 
Sint 

TeOSPCoswy 

O cos¢sind 
Sinw 

—Singdcosy 


S08 


the total trans— 


denotes its 


a 35 =—SannG 
Gda-.= Sindcose 
a 
(8 ) 
on cOsgcosdO 


In relating the angular velocity rates of the aircraft, 


p, gq and r to the inertial frame, we make use of the trans— 


formations one eyo. and 


Va 


This development 


Steals projecting the orientation angular rates, w, 6 and 


Memento the aircraft body axes. 


In addition, the inertial 


angular rate (y) must be included EO e2CcOumE stor “Ehe -uOta— 


Bren Of the carteSian coordinate reference frame (X,Y 


Both y and i ase considered tombe: applied about the Z 


page 


ee 


1 


axis of the cartesian coordinate reference system and, 


eicough th 7 
: FN p7 i 


mealized. 


volves only two transformations 
jection on the body axes. 


the body axis involves no transformations. 





expressed: 
Pp 0 
rele 0 
v U+y 


} transformation, 


The Euler pitch rate (hme aboueetne . 


IT 372 


}) and 


their projections are 


5 axis, in- 


IT /3)° 


FOUN Dro 


The Buller roll rate (6) about 


We 


Mathematically 
> 
+ 0 fo) 
0 





Marsubstitution of Equations (5), (6) and (8) for [T3 72). 


and [T respectively, in Equation (9), the exact 
IT 73! yal p vi q 62) 
relation between the aircraft angular rates and the Euler 
angle rates is given by: 
p= ¢ — sind (W+y) 


6 cos¢d + Bee ebiie Caay (ee@ ) 


q 


i cosécosé (Wty) —§ Sind 


The time rates of change of the direction cosines Se 
and the angular rates (p,q and r) are obtained by direct 
differentiation and are presented in Appendix A. 

Pee sAGRANGE FORMULATION OF THE AIRCRAFT EQUATIONS 

OF MOTION 

Assuming that the aircraft 1S a Single rigid body, the 
three position coordinates (ZorRrY) and the three orienta- 
tion angles (W,8,¢) can completely describe the six degrees 
of freedom of the vehicle. These Six variables form the 
independent set of generalized coordinates paeeeeee fOr 
Lagrange modeling of the motion. Let T denote the kinetic 
energy of the system; V, the system potential energy; q_, 


ie 


megeneralized coordinate; and F a generalized force due 


dy 
to non-conservative forces. Then the Lagrange equations 
for this system can be written: 
so) - St Ba rg, v= 1,2,.-06 (11) 
oq dy dy 


a 
the Kinetic energy of the aircraft, relative to the 


inertial reference frame, is given by the superposition of 
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the energy of the center of gravity, plus inertial rota-— 


meena terms about the center of gravity, and is written: 


Meee. 2 «2 = 2 7 2 2 2 
= aml (YR) ve el) A (20) lla er 5 1Iyp ee ciaa Ee }+ vos (1 2) 


It is necessary to use the relations given in Equation (10), 
as EGQuation (12) is not compatible with the generalized co- 
ordinates (ZeRs¥ W079). Upon substitution of these rela-— 


tions for p, gq and r the compatible kinetic energy equation 


nS: 
T = sm{(yR)*+ (R)*+ (2.)°) + ZI, 1b — sine (j+y) 17 
+ sty {6 cos¢ + sinécosé (Wty) ]~ (203) 
+ =I, [cospcose (b+y) 1° 
+ eo —siné (W+y)] [ (Wty) cos¢cosé — 6 sind] 


If a reference for the potential energy is established 


at sea level, V 1s given by: 
V = mgZ (14) 


The partial derivatives of Equation (13) with respect to 


the time rates of the generalized coordinates ga ) are: 
aq 
oT . 
92 : 
O 
=) mk (16) 
dR 
pe mRyt+ icra. | ga,+ Lara 
: KS Vo 23 L335 


oY 


+ (op cosdcos6 — r sin®) (17) 


Iyg 
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and 


aT 
ay 
ar 
36 
or 
36 


The 


~ TPO, 3 oc T9053 - Ira, = TVoP Sind (18) 


r Sino + I 


I! 


Iq Cosa + I, yz (pcos¢coso—rsind) (19) 


T.P ae (20) 


xz4 _ 


time rates of change of Equations (15) through (20) 


represent a qa ) and are: 


dt 
dq. 
d aul s 
a (21) 
Gt 5% O 
O 
- (27) = mR (22) 
dR 


= MR(Ry+2Ry) +1, (pa, 3tpa,,) tl, (qa,,+Ga,3) f 


iE (ro 


z, tr.) +1,, (—p(68cosdsin6+¢singcos6) + 


53 


~ 


peosdcos6— @rcosd — r sind) (23) 


= I, (pa, ,tpa, 4) +1, (qa,,t+ga,,) +I, (ra,,+ra,,)+ 


I, (-p (6cosdsind+psingcos6) +peospcosé = 


6 r cosé ~-r samo) (24) 


aes ce I, (qcos¢—-q$sin$) +1, (rsing+rgcoss) 


d_ (at, 
Cte zy 
Gea,.o'T 
ae) = 
Qt a0 
d ,3T 

068 
Gee nol 
Beene (—) re 
at a6 


—I,, (pocosot+psing) 20) 


= Typ 7 1,54 (26) 


ae 





The partial derivatives of Equation (13) with respect 


to the generalized coordinates ee ) are: 
ie 
oT 
O 
oT 2 
<= = my°R (28) 
ak 
a. 
a =n C30) 
ah ee 
ae (p+y) (IT ,pcosotiygsindsingd+I rsindvcos®) 
+1, (—pcososiné (pty) a rcosé (Wty) ) . (eon) 
a = qr (1,-I,) + I, (-cos¢ a sindcosé (Wty) ) (3:2) 


The partial derivatives of Equation (14) with respect 


to the generalized coordinates oo ) are: 


1a 
ie , 
O 
av _ av _ av _ av _ av _ 
) ; Snr, a 0) > ‘ ey 


Memataens (21) through (34) complete the formulation 
of the left hand sides of the six Lagrange equations, as 
generally expressed by Equation (11). In completing these 
equations it is only necessary to formulate the generalized 


forces FL, FP F By Penida es Thnesorinciple Of yvir— 


a Re y wy’ Qj ¢ 
miciework forms the basis for these calculations and can be 


described mathematically by: 
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J€ 


OX oY X 
oq. ) + Py ‘9q_ ) + Fy (22 ) + eG.” +4 
dy iS 
JE, JE, 
Ny (Sq ) + Mo(a—") = Fa, (35) 
ee : 


where X,Y,Z are positive incremental displacements along 


the X axeS, € ,€.,€ are positive incremental angu- 


Me?) x EyrEe 
lar displacements about the X11¥4 12, 


The incremental translation or rotation along or about 


ase Sin 


ae or ay axis, caused by varying each generalized 


coordinate individually, can easily be visualized with the 


the X 


aid of Figure 3. For example, varying the generalized co- 
ordinate Z, by a positive increment will produce no trans— 


1 Or Yue ana no rotation 


about Z,- WeormmowEgGuiatton (35) as a guide (realizing that 


each X,Y and Z appearing should be subscripted with a one, 


maeron Or rotation about either X 


because we are Ron with the carteSian coordinate refer- 
IX, oY Jey Jey, 
ero frame), it follows that a7. ' RZ! <5 ; wa 
€ O O O O 
57 ~=OCOarre idemteloatily Zero. From Figure. 5 1 1S clear that 
O 


and 





Figure 5. Positive Incremental Displacement of the 
Generalized Coordinate Zoe 
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as Z is incremented a positive amount, 2, follows direct- 


al 
fieeeeits megative defined direction. This implies that 
dZ | 
— —l and the generalized force Fy = —F, follows from 
O O ¥ 


Mewation (35). 


The remaining generalized forces with respect to the 


cartesian coordinate reference frame (X1,Y 74,) can be 


derived in a similar manner. The complete set is: 
F, =M 30 
p> Ms, (36) 
Fo = —M, Sin wy + M, cos w (3:7) 
1 a 
ae = My cos0cosy + My, cos@sinw — M, sind (38) 
il i il 
> -F, (39) | 
O 1 
R Ya 
foe FF. R+ M (41) 
Y Xy An 


Making use of the axes transformation Eey ale? the general- 


ized forces expressed in terms of the body axes are: 


- = (Oj) oF y + Aa5F, + A25F,) (42) 
Es = “R(a))F,+ SM legs Phe, Ral UCR satel i eee (43) 
Fy = 013M, + 453M, + a,5M, (44) 
wee fh, cos } — Me sin 9 (45) 
Me = My (46) 
a em a oat y * Yah) se 
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Adding Equations (42) through (47) to the respective 
meant —nand Sides of the earlier mentioned grouping (of the 
left-hand terms), in the general expression given by Equa— 
tion (11), completes the development. This grouping is 
shown in Appendix A, where the six equations of motion are 
presented in their entirety. 

Engine thrust terms are not included in the six equa- 
tions of motion. Previous studies indicate that use of 
thrust in the spin is normally avoided since flameouts are 


likely and serious damage can result (Ref. 4). 


PeeothEADY SPIN EQUATIONS 

The full equations of motion (as given in Appendix A) 
can be greatly simplified for a steady state spin condition. 
mimcuch a State, the aircraft 1S imagined to be descending 
at a constant velocity, at a constant radius (R) and rota- 
meon rate (y) about Ehemcentral Spin axis and with a fixed 
errentation with respect to the carteSian coordinate refer-— 
ence frame (Xj) /Y,72))- This implies that all second-—deriv- 
ative terms and at first-derivative terms of the Euler 
angles (W,8,¢) can be set equal to enon On examination 
of the auxiliary terms of Appendix A, it is easily seen that 
meet ime rates of the direction cosines an and angular 
accelerations (p, q, r) are equal to zero. After making 
Bese Simplifications, the six equations of motion that de— 


Scribe the steady spin condition are: 


pe (1) 2 0 043Fy, + GooFy, + 45,F, — mg (48) 
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RES (2) = 0 = Qs oF, + AooFy + A.5F, — my R (49) 
RES(3) = 0= R(a,7 Fy to5,Fy + a5F,) + 
(453M, FO5 My + a.5M,) (50) 
RES (4) = 0 =a, 3My + a5 3M, ua a.3M,) (5) 
Pests) = 0 = —y (I,pcosé + I,gsingsin6é + 
I, (rcos$sin@) + M,cos$ _ M, sind (52) 
-y1,, (pcososind + rcos6) 
RES (6) = 0 = My-ar(I,—I,)—-I,,pysing¢coso (53) 


mnie aerodynamic coefficients are written in a form con— 
sistent with the data obtained from Ref. 5. The aerodynam- 


moemuorces and moments are written: 


oo al 2 e = 
ge SC, + Cx, Spc Cx q) (54) 
2 S = = 
Fy = 50V pry y, ay 6 tC, ‘p+ ge (55) 
ine 2 ie r = 
ee oY S(C, + - on sans cl) (56) 
i e = Z 
My = 50V Sb(C,+ Cy ©6 + ee ae, ap (a7) 
- 6 Pp Yr 
ie — ‘. 7 58 
My = ZV Sc(C_+ Cm « oe + “mV (58) 
e = — 
M_ = i V-SbI(C +C °*6 +C TO. G p+. er) bso) 
L 2 ie ee) ere Tei n 
Oo. 6. p 08 


iPhmthe relations given by Equations (54) through (59) 
were substituted into the steady spin equations, as given 
by Equations (48) through (53), the six equations would 
then describe the resulting aircraft state in terms of con— 


ventional aerodynamic forces and moments. 


Zor 





The residual notation (RES(i)) introduced in the steady 
Spin equations, as given by Equations (48) through (53), is 
used aS a measure Of how closely each associated equation 
is being satisfied. Ideally, each residual should be iden— 
tically equal to zero but this is not possible due to trun- 
Seren and round—-off error within the computer at various 
stages of tne program. 

Peecriteria function 1s defined that will indicate how 
close the complete set of Six equations is to a steady spin 
Geme@qation. The criteria function chosen for this investi- 
gation was the sum of the squares of the individual resid- 
uals, which is written: 


F = RES(1)°+ RES(2)°+ RES(3)7+ RES(4)°+ RES(5)° + RES(6)~ 


(60) 
meeseruncckion will be the object of the minimization effort 
in the optimization procedure described in the next section. 

The airplane has six degrees of freedom with second 
order differential equations; therefore, 12 variables are 
required to eeceree its state. For convenience the 12 
variables chosen are: VV 106,056,252 ,R,R, y and Y. Mak—- 
ing use of the fact that each equality constraint imposed 
upon the system will remove one variable, the total number 
of required variables can therefore be reduced. Some of 
these have been discussed in describing the steady spin; 
Specifically, that R = w = § = b = 0. They were also used 
to reduce the general equations given in Appendix A to the 


steady spin equations given by Equations (48) through (53). 
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By specifying that the search be made at y = 0 and at a 
specified altitude (2), the total number of equality con— 
straints 1S six. This leaves Six variables to search over 
meweach control surface is held in a fixed position. The 
chosen set of variables to be searched over in this inves-— 
Eegation are: Y, Re Zor Ue amend o. Specifying these six 
impetal conditions, the altitude parameter (p, as Zo 
doesn't appear explicitly in the steady spin equations), 
the phySical parameters of the aircraft (mass, chord, Span, 
etc.), and having a set of tabulated aerodynamic coeffi- 
ements, the criteria function can be evaluated. MThe only 
thing now required is a systematic method of varying the 
S1xX variables (yY, R, Zo We.) tO ObEain, Che minimum 
meme OL the criteria function close to zero. 

The solution of the steady spin condition involves the 
Simultaneous solution SiMsisenignly coupled nonlinear dif— 
ferential equations involving aerodynamic data. There 
does not exist a closed form analytic solution to this set 
of equations and therefore a numerical method must be 


employed. 
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meer VvlotD FORMULATION OF THE STEADY SPIN EQUATIONS 


The steady state equations were FORTRAN coded for use 
mmecOogrameSPIN and anitial efforts at predicting the equi- 
librium spin condition were unsuccessful. During these 
efforts it was realized that the rate of convergence was 
extremely slow and that the solutions were not close to the 
states already known to exist, as shown by Adams (Ref. 7). 

After some analysis it was theorized that the deletion 
of the acceleration terms from the derived equations of 
motion was the cause. The equations should have been un- 
coupled and solved explicitly for these acceleration terms. 
Then the residual terms, representing the acceleration 
values, would be a true measure of just how close the air- 
See 15 to the equilibrium spin condition. 

Attempts at sheizonre eae thestulls equations Gf MOtLOn, as 
given in Appendix A, showed that it would be, at most, a 
very time consuming and nearly impossible task as the yj, @ 
and ¢@ equations are so highly coupled. Instead, a search 
for an alternate set of moment equations was made with the 
decision to use the equations given by Etkin (Ref. 3). The 
three equations express the moments about the body axes and 


ate 


M, = Iyp — oe rs gd — I,,Pq (61) 
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7 : _ a g 
My = cincl “p rp(Iy + Io (Pp q) (62) 


M =-I, oP + I,v + pq (I — I,) + J (63) 


Z, xzor 


The first three equations of the original formulation 
express Newton's second law of motion along the radial, 
tangential and vertical directions through the aircraft's 
center of gravity. These equations can be expressed, 
without introducing the transformations neve about the 


carteSian coordinate reference frame as: 


F. = my (64) 
al 
oe 
F + my R = mR (65) 
Y 
fl 
F = mz (66) 
aa 


These three force equations when combined with the three 
moment equations, as given by Equations (61) through (63), 
describe the motion of and rotation about the center of 
gravity of the aircraft and therefore completely describe 
the motion of the six—degree—of—freedom vehicle. 

Bye explicitly solving Equation (64) for the accelera— 
mon term Y, Eola clOn G5) ) for R, Eouatton, (60) “for Z and 
Pemation (61) through (63) for D, gq and r the equations 


are in the form desired for the optimization search. 


a 





ey Or | IMP aAT TON PROCEDURE 


CONCEPT OF OPTIMIZATION 


Using state variable notation, the equilibrium solution 


and optimization scheme for this particular formulation can 





be described mathematically as follows. If we let 
R 
Z 
* | = ; , the acceleration terms 
p 
q 
r 
and 
y 
it 
Z 
1: | = , the independent variables, 
y 
6 
p 
therefore | s | = ) Py . Then the optimization scheme 
mind S ! x | euch that x ! =O GOr «equilib lis The Op 





timization scheme employed in this investigation uses the 


criteria function which is the sum of the squares of the 


acceleration terms (modified equation residuals), given by: 


J(¥) = 


ee ee Re leas 4 ee 
i=1 


a2 





meaca Linds | xe such that GY) is an absolute minimum 





= 0. 


which yields ! x 





B. OPTIMIZATION ALGORITHM 

The optimization method used in this investigation 
menes the Minimum value of the criteria function, which is 
a real-valued function of several Teen Ther criteria 
Mumetion is used as an indicator of how close the aircraft 
mS)tO being in an equilibrium spin condition. The minimi- 
zation process takes into account any constraints imposed 
upon the system, which, in this case, were imposed limits 
of angle of attack and sideslip from the available aero- 
dynamic coefficients. 

The optimization procedure starts with a given set of 
independent control parameters and varies the current values 
of these parameters, within the imposed boundaries, to im—- 
Breve the value of the criteria function. The process of 
parameter adjustment continues until the minimum of the 
@eeeeria function 1s reached. Depending on the hyper-sur— 
face defined by the independent control parameters, the 
Minimum found may be just a local minimum which would not 
Satisfy the slesliaeel steady—-spin criteria functional value 
(equal to zero). The absolute minimum of the hyper-surface 
will yield a criteria function value equal to zero and is 
Paewonyect Of the search. The initial conditions given 
for the independent control parameters will greatly influ- 
Peeesthe Success of finding the absolute minimum. There— 


Bere, there is a certain trial and error process in the 
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search for the absolute minimum without any guarantee that 
one exists. This then is the main disadvantage of the op-— 
timization procedure because in spite of the energy and 
Memeoecverance Of the researcher he cannot, at the end, say 
that a steady spin condition does not exist for that air-—- 
craft if in fact he hasn't been able to locate an absolute 
minimum. .This 1S a very common disadvantage shared by all 
optimization procedures available at this time. 

aie Optimization algorithm used in the spin program is 
called EXTREM. A complete listing of this subroutine is 
given in Appendix B along with a functional flow chart. The 
user specifies a set of independent variables in name and 
weeralt value and a corresponding number of scaled incre— 
ments defining the initial step sizes along each of the in- 
G@ependent variables. The algorithm first checks that the 
Seven initial conditions do not violate any constraints 
imposed on the system and then proceeds to determine the 
main line of search in the variable space. The initial con- 
ditions form a point in this variable space and the speci- 
fied initial step sizes give a second. A vector joining 
these two points determines the initial main line of search 
Perecne first stage, as shown in Figure 6 for a two-variable 
problem. Along this main line the approximate minimum point 
1s found by extrapolating a parabola through three points 
- DX, X,. and X.+ DX and then interpolating 


a 1 


“A Gram—Schmidt orthogonalization. 


determined by X, 


fot the Minimum (X44). 


Process determines the secondary direction of the first 
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stage which is orthogonal to the main direction at Xo 44° 


The approximate minimum point Xs 5 to eround along this Line. 


SECONDARY DIRECION 


(2) OF FARST STAGE 








y Maw DIRECTION 


\ 
a AY OF FIRST STEQE 
Xray dd 


A So 
—o aN 
Be Gece, Aqaz 


j 
(GUESSED INITAL Pont ) 


X(!) 
megure 6. Main Line of Search in a Two Variable Problem 


Meme third variable (not shown in Figure 6) a line of 
search through ee is determined by the Gram-Schmidt pro- 
cess that is orthogonal to both the main and secondary 
directions. Each new direction is orthogonal Ome ela the 
preceding directions of search and the procedure for deter-— 
mining these lines of search and approximate minimum points 
meecontinued until the minimum along the ae line (K being 
the number of independent variables) is POUMG area core mass 
point the first stage is complete and the main direction 
for the second stage is determined by a line joining the 
point Xx, and re where Xia is the approximate minimum 
along the ee line of the first stage. The second stage 
proceeds just as described for the first, and the overall 
stage procedure continues until a local or absolute rela- 


tive minimum is found. Between stages the program checks 


aD 





mae stOpping conditions on the criteria function value, 
the arguments, and the maximum number of stages allowed to 
reach the minimum, and will stop if any of these conditions 
is met. The step size along a line of Search is variable 
in the program, which allows the search steps to decrease 
near the minimum or within tight curves and to increase 
again once these curves are past. 

Peemore AGetalled description of the optimization method, 


the algorithm and several test cases can be found in Ref. 8. 
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V. DISCUSSION OF RESULTS 


fee VERIFICATION EFFORT 

Iittean attempt to verify that the spin equations of 
motion describe the Ore Morton and also to check out 
the entire computer program, several runs were made for an 
equilibrium flight condition of straight and level flight. 
The choice of this flight mode was made so that the results 
could be verified by hand while searching over a reduced 
number of independent variables in satisfying the same cri- 
teria function as discussed earlier. 

As mentioned earlier, thrust terms were not included in 
the general equations of motion, aS previous studies indi- 
cated that its use 1S normally avoided, since flameouts are 
likely and serious engine damage can result. However, in 
order to simulate straight and level flight thrust was in- 
cluded with the very restrictive assumption that it was 
Orientated in enaunecaeive X direction and passed through 
moe center Of gravity. 

Straight and level flight can be simulated in program 
SPIN by specifying a value of equal to 90°, y and Z both 
Set to zero. By assuming that a wings-level, straight line- 
Seer light Orlentation can be obtained with neutral rudder 
and aileron control, the value of 6. and 6 are set equal 


to zero as well as the bank angle ». The value of the 
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mieinds Vector (R) 1S arbitrary. With this accomplished for 
each run a value of the pitch angle (8) (in this special 
case @ is the same as the angle of attack) 1s specified. 
Three independent variables (thrust required, elevator 

angle and airspeed) and their initial step sizes are speci- 
fied in the main program ceic routine yields the value 

of these three variables that would be required for straight 
and level equilibrium flight. 

The results of these runs were verified by uSing an 
electronic calculator; the equations of motion were satis-— 
fied. The numerical results of thrust, velocity and eleva- 
tor deflection angle were reasonable for that particular 


mmecght Configuration. 


Peer ROGRAM SPIN RESULTS 

Several SPIN runs were Berueeed or Conte igubration. 8 in 
order to compare the results with those given by Adams 
(Ref. 9). The mass, dimensional characteristics and aero-— 
Semamic coefficients for this aircraft configuration are 
listed in Appendix C. 

The results of the predicted spin characteristics for 
configuration B are shown in Table I where they are com— 
pared with the results obtained by Adams. The predicted 
spin characteristics developed by program SPIN, using the 
optimization algorithm EXTREM, are in good agreement with 
those predicted by Adams. The slight difference in the 
Meese, second or third decimal digit in some cases is prob— 


ably a result of the particular computer being used. 
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Adams used a CDC machine with a bigger bit word than the 
IBM-360 and therefore more significant figures could be 
handled in the numerous calculations performed. AS a re- 
pmeewot this the minimum value of the criteria function 


obtained with program SPIN was Mone. Adams has been able 


memorive his criteria Benen on residuals (also a sum of 
time squares but of a different set of terms) to a value of 
moe ° 
Mie third grouping of results in Table I show very 

little similarity, despite the same control settings being 
used. The first steady spin predicted is an entirely dif- 
ferent equilibrium condition as can be seen by comparing 
the values of the radius (R) and angular rate (y). The 
second spin predicted would at first appear to be the same 
as that given by Adams. The difference in values of the 
angle of attack (a), inertial velocity (V) and radius (R) 
make it an Bats GQiiferene Soin CcOnadlt1on and not 7ust 
a poor representation of the same spin. This statement is 
justified by the pilocen correlation of results, in the 
other groupings, for the same value of the criteria func 
mien (10 -°). | 

ie COmputer prediction runs for configuration A are 
Shown in Table II. Ee eee no other results to compare 
with this configuration. The spin prediction scheme was 


meplied to this particular aircraft to assure that the 


Temidod would at least work for another aircraft. The choice 


oo 





FOr this second configuration, over other possibilities, 
was made easy because the aerodynamic derivatives were al- 
meady in tabulated form as given in Ref. 7. 

The effects of density variation on the equilibrium 
Spin conditions is shown in Table III. Grantham and Graf- 
mom (Ref. 8) concluded eereenieter Studies that a decrease 
in density gave a faster rotation rate (y), higher veloc— 
ities and little change in angle of attack and sideslip. 
The parameters predicted by program SPIN verify this be- 
havior with the added effect that the radius of the spin 
tends to decrease, 1.e., the spin becomes tighter with a 
decrease in density. The change in the orientation angle 


86 would indicate that the spin would be somewhat flatter 


with a decrease in density. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


As a result of developing and presenting the analytic 


spin prediction technique, the following conclusions and 


moeservatlons are made. 


Ie 


_— —-— 


The original developed equations of motion, as 
presented in Appendix A, are too highly coupled 
memeemal Crrecrtive tool in this type of prediction 
scheme. In this light, those equations can be con- 
Sidered, at most, as an academic exercise leading 
momomerather uUnlgue formulation. 


The comparison of results demonstrated the accu- 
racy of the method and the effects of density 
variation study showed one possible use of the pro- 
gram. The capability of this method to predict 
Steady spin conditions is limited only by the accu- 
racy with which the mass, phySical characteristics 
and aerodynamic coefficients are determined. 


Buer eo une dlsece modeling in a cylindrical coordi=— 
nate reference frame this method has a certain 
advantage over that demonstrated by Adams (Ref. 7). 
The number of equality constraints on the solution 
required by this prediction method is far fewer, 
which allows for a much simpler computer formula-— 
tion and equally accurate results. However, this 
method uses more computer time than Adams' in 
ConvengingetOoewene equilibrium Spin condition. 


As with any optimization technique the importance 

of equally weighing the residuals of the equations 
is of utmost importance. The value of the residuals 
should be as close as possible to each other, in 
order Of magnitude, for proper convergence of the 
optimization algorithm. Unequal weighting of any 

of these residuals cause these terms to dominate 

and influence the search in a direction that doesn't 
actually contain the desired solution minimum. Many 
hours of research time were lost by not recognizing 
thes InDpOrtant Lact . 


In light of the demonstrated usefulness of the analytic 


prediction scheme the following follow-on work could prove 
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valuable in better understanding the spin problem. It is 


recommended that: 


Ie 


The use of rotary~-balance aerodynamic data be in- 
vestigated as a possible way of more closely relat— 
ing the aerodynamic coefficients of the tunnel 
model to that experienced by the real airplane. 

The use of this type of data will allow at-—desk 
CalleculaeloneOrerne Steady Spin conditions. This 
particular spin estimation method, if used to pre- 
MiettweniesIMledal COnditions for program SPIN, 

would reduce any senSitivity to the guessed initial 
conditions and also reduce the required computer 
convergence time. This particular method of spin 
estimation is given in Ref. 9. 


The equations of motion could be linearized in 
order to obtain some information about the stability 
CMmeneOmearGted Equilibrium Soin condition. 


This prediction scheme could be coupled with exist- 
ing fixed-base spin simulators to provide the 
Starting initial conditions for the spin simulation. 
Wiehe tiio the recovery characteristics could be 
investigated, recovery techniques perfected and the 
necessary pilot interactions evaluated and studied. 
This one tool alone could provide valuable insight 
into the problem areas and shows much promise in 
future work. 
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PePEN DEA A 
PULLS BOUATIEONS OF MOTION 


The following six egquations completely describe the 


bigid body motions of an aircraft, using body axes, as 


formulated in a cylindrical coordinate reference frame. 


Z EQUATION 
ma aa 053F%y as Q53F. ts A53F, + mg 
R EQUATION 
oe 2 
Mee my R = (a, 5Fy + A55Fy + a55F,) 


y EQUATION 
mR(Ry + 2Ry) + Ty (po, 3 a po, 3) + I, (qa,3 + Ga,3) 


+ I. {-—p(@cososind + osindcossé) + 


AG ‘ 
p cos¢cosé’— 6rcosé — rsind} = 
ke iera cone 2 ag? 2)) Saal P O5aMy tO 35M 


» EQUATION 
I, (pa, 3 : pa, 3) + I, (qa, + qa, 3) + 


I, (ra,, + YO. 3) + IL, i-p (8cos¢sing ~ 


osingcosé) + p cosgcosé — 6r cosé — r siné} = 


ee oa wy) aah 
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9 EQUATION 


i (g COSda.- qosind) _ I, (rsing + rdcos$) 


+ (y+y) {I pcosé + Igqsingsin® + I,rcos¢siné} 


—I,, (¢cos>p apeeseleT) DO.) ert TL, (pcospsind (pty) 


+ rcos0 (W+y) ) = Mee COS. M, Sing 


xy 

@ EQUATION 

TP t q, (1, = I.) _ I, (q + p(¢cosd + sindgdcosd (pty) )) 
where, 

1, = cosycosé@ 

or = —Wsinycos6 = Acosysind 

Ai5 = Sinycos9§ 

cae = Woosycosé = Osinvsiné 

53 = —=Si ng 

053 = =OCOSG 

Qo, = cosysinésingd — sinycosd 

oe = —Yysinvsindsing _ 6cosycosésind 

—bcoswsindcosy = cosycosé + osinwsing 


Sinysindsingd+ cosycosd 


woosysindsind— ésinycosésingd 


— osinysinécosd — wsinycos¢ a scoswsing 
cosésind 
pcosécosd — ésinésing 


cosysin@cosd + sinysind 
—isinvsinécos¢ + Bcosycosécosd¢ 


—odcossin@sind + wooswsind + ASinYCcosd 
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32 


53 


a3 


pe ae 


‘(om 


SinWsin8cos¢? — cosipsind 


WooswsinOcosw + @sinycosdcosé — psingdsinOsing 
+ Wsinysing = coscosd 

cosdédcos® 

= Asindcosd a bsindcosd 


b — sind (W+y) 

$ ae wOcosé = Wsind = yOcosé = ysino 

Acoso + cosésind (Wty) 

= 6osing + 8cos¢ a wOsindsing = WdocosdOcosé 
+ oos@singd — yésingsind + ydcosOcos¢ 


+ ycos@sind 


cosé cos¢ (Wty) — @sing 
— wAsinécosd — wocosdsing os wWoosdcos¢ 
rt Adcos¢ = Osing = yésindcos¢ 


— yocosésing + ycosécosé 
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APPENDIA.B 
PROGRAMS SPIN 


A. INTRODUCTION 

The intent of this eon 1S to’ serve aS a users 
manual or guide to the use of Program SPIN. An alphabeti- 
cal listing of the computer variables and their meanings 
1S provided. A brief description of the main purpose of 
each subprogram aS well aS an extensive explanation on pre- 
Pamang the data deck is included. A listing of the program 
with numerous comments is included along with several func— 
tional flow charts to help in understanding the structure 
and relationships of the subroutines. A sample output is 


given and the storage and time requirements of the program 


are discussed... 


Eee olSt OF COMPUTER VARIABLES 
The following. 1S an alphabetical listing of the compu- 


ter variables and their meaning. 


Er l=A 33 Direction cosines 

ax) Working matrix for Subroutine EXTREM 
ALPHAD Aircraft angle of attack, degrees 
ALPHAR Aircraft angle of attack, radians 

B Wing span, ft 

BETAD | Angle of sideslip, degrees 

BETAR Angle of sideslip, radians 
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CREAR 
CL 
CM 
CN 


COEFA( ) 


SOrnrAB( ) 


CPi 
CPSI 
CR( ) 
CTHETA 
CX 

CY 

CZ 

D 


DDELA 


DDELE 


DDELR 


DELTAA 


DELTAE 


DELTAR 


DF MAX 


DGAMD 


Wingechora,. ft 

Otter moment coefficient 
ToOtatawertch moment coefficient 
Total yaw moment coefficient 


Aerodynamic coefficients which are a function 
Gio? LeomOnrmcdetack. Ol Vv 


Aerodynamic coefficients which are a function 
of angle of attack and angle of sideslip 


Cosine of PHI 

Cosine of PSI 

Working vector of interpolated coefficients 
Cosine of THETA 

Total X force ecoerricient 

Total Y force coefficient 

EOtal LZereonece COCELICIEnNE 


Oma) 27 KE * Zo 


Incremental step size of aileron control, de- 


grees 


Incremental step size of elevator control, 
degrees 


incremental Step Size of rudder control, degrees 


Aileron control deflection, degrees (Positive 
with right aileron trailing edge down) 


Elevator control deflection, degrees (Positive 
with trailing edge down) 


Rudder control deflection, degrees (Positive 
with trailing edge left when viewed from above) 


SEOpeINGmecoOnadltl1on Of the functional variation 
in Subroutine EXTREM | 


Incremental step size of GAMDOT, rad/sec 


a0 





pen rT 
wes I 
DR 
DRDOT 
BaHETA 
DTHRUS 
DXMAX 
ee ) 
DZDOT 
is 


FMAX 


Eee | 


GAMDOT 


HRHOS 


IW 


LMAX 


NERINT 


PBAR 


Incremental step size of PHI, degrees 
Incremental step size of PSI, degrees 
Incremental step size of R, ft 

Incremental step size of RDOT, ft/sec 
Incremental step size of THETA, degrees 
Incremental ance Ore TA RUST EDs 
GOnesesoarametern For Subroutine EXTREM 

Array of step sizes of the indeecncens variables 
itiememental Step size of ZDOT, ft/sec 


HuneeLon to be minimizea 


Desired minimum value of the optimization 
fume & LON 


Optimum value of the function being minimized 


CamminatLenakl acceleration, ft/secz 


Medteomonuewange Of the cylindrical orientation 


variable y, rads/sec 
0.5*RHO*S 


COuer Ole Darameter OF SUbrOULELNe EXTREM for 
Weenie struet Ons 

+1 All output suppressed 

Peeinaly output Oly 

tomeOUcoUE at tne end Of cach stage 


The number of independent variables being 
searched over 


SeOpulnc COndleLlon on the number of stages Of 
Subroutine EXTREM (A negative sign indicates 
that a minimum is sought) 

Print indicator, set in MAIN and passed to 
Subroutine DATAIN to echo check the aerodynamic 
coefficients ( =0, No print; =l, Prints) 


Angular velocity about the X body axis, rad/sec 


PAB /2V 


ok 








eis 
BAIR 
Pol 


ow R 


QBAR 


RBAR 


reo Ty 


RES ( 
RHO 
RZ 
S 
Seit 
Seo lL 


SSAM 


STHETA 


SUMRES 


THETA 


THETAR 


V 
VSQ 


VXB 


VYB 


Euler orientation angle ¢, degrees 
HibeCrerGnrentdtaon angle 46, radians 
BUlerENOntentation angle wy, degrees 
Euler orientation angle Ww, radians 
Angular velocity about the Y body axis, rad/sec. 
Q*CBAR/2*V 
Chine 1 Cal@eeadanis seer t 
Ra 2 ov 


Rate of change of the cylindrical radius, 
EE /SeC . 


Residual of any of the spin equations 

Density OF alr, lbs-sec*/£t‘ 

Angular velocity about the Z body axis, rad/sec 
Wing area, ft- 

Sime Wot sy PHI 

Sine of PSI 

Aircraft mass, Slugs 

Sine of THETA 


Sum of the squares of the individual equation 
residuals 


Euler orientation angle 8, degrees 

Euler orientation angle 6, radians 

Velocity of the vehicle center of mass, ft/sec 
Velocity Sauaen 


Vehicle velocity component along the X body 
axils, ft/sec 


Vehicle velocity component along the Y body 
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XFORCE 


XI 


XMOM 


nORCE 
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YMOM 


ZDOT 


ZEORCE 


a1 


ZMOM 


2X1 


VYB/V 


Vehicle velocity component along the Z@ body 
axis, ft/sec 


VZB/VXB 


Working vector that contains the independent 
variables 


Aerodynamic force along the X body axis, lbs 


Moment of inertia about the Z bedye axis, 
slugs-ft? 


Aerodynamic moment about the X body axis, 
ifPe— Wis 


Aerodynamic force along the Y body axis, lbs 


Moment of inertia about the Y body axis, 
slugs-ft2 


Aerodynamic moment about the Y body axis, 
ieee 


Vebuele velocityain the cylindrical axis 
dmrectlon, —tescec 


Aerodynamic force along the Z body axis, lbs 


Moment of inertia about the 2 body axis, 


slugs-ft? 


Aerodynamic moment about the Z body axis, 
Et=ilbs 


Cress pProauct Of inertia about the xX and Z 
body axes, slugs-ft¢ 


See PROGRAM DESCRIPTION 


Peog@ate oft CONSIStS Or ja main controlling program and 


Six subprograms as illustrated schematically in Figure Bl. 


hice Mir pregran S, most important function 1S to set 


the independent search variables and their initial step size. 


Otherwise the MAIN program is basically an input/output 


aS 





mieine. The initial conditions and optimization control 
Barameters are read into common memory. The initial con- 
@re1ons are printed out aS are various final conditions 
ee the airplane. 

Subroutine DATAIN is called by the MAIN program to read 
in the aerodynamic Pee oon Coefficients that are a 
jmometion of both alpha and beta are read into the COEFAB 
mieay and those that are a function of Piene only are read 
mee the COEFA array. DATAIN will write out all of the 
MeerticientsS if the variable NPRINT is set equal to one. 

Subroutine EXTREM is the optimization program which 
determines the minimum of the criteria function without 
levying tO Calculate any derivatives. The single most im-— 
portant function of this subroutine is in determining the 
optimum direction of search since the criteria function 
being minimized resides, and is evaluated, in Subroutine 
FN. The returned values of the criteria function are used 
feminiterpolate a parabola along a line of search, to 
determine the ae eeu on this parabola and to determine 
new directions of search. 

enero tine FN evaluates the criteria function value 
=emGeough Calls to Subroutine CONST, Subroutine CLOKUP and 
Subroutine SSEONS. In addition, boundaries on the search 
are -imposed due to’ the limited range of angle of attack 
and sideslip angle aerodynamic coefficients available. 

Subroutine CONST calculates various equation parameter 


values needed in the equations of motion. Values of 
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angle of attack and sideslip are determined so that the 
limits of the available coefficients aren't exceeded. 

Subroutine CLOKUP takes the passed values of angle of 
attack and sideslip and performs a linear interpolation 
on all the aerodynamic coefficients. 

Subroutine SSEQNS Seneanne the six equations of motion 
of the airplane equated to the residual (RES(I), I=1,6) 
vector. The residuals are evaluated and passed, through 
common memory, to Subroutine FN so that the criteria func— 
tion can be evaluated. 


Functional flow charts of the MAIN program and each 


Subroutine are shown in Figure B2 through B8. 
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SUBROUTINE 
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Hg uress 1. 


MAIN 
PROGRAM 


PRIN 
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SUBROUTINE 


EXTREM 





SUBROUTINE 
FN 
SUBROUTINE SUBROUTINE 
CLOKUP SSEQNS 





Schematic of Program SPIN 
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See 
Peet Sew omARkeH VARIA BIIE SEEP INCREMENTS 


READ INITIAL CONDITIONS AND 
PATREM CONTROL PARAMETERS 


Se eC Onr sere Ni Phin IND CATOR 
Chl the NOMBEREAND SPECIFY THE 
SEARCH VARIABLES 






CALE SUBROUTINE DATALIN 


PRINT THE INITIAL 
CONDITIONS 


CALL SUBROUTINE EXTREM 








PRENT The puna 
PARAMETERS AND 
RESIDUALS 


| STOP 


Figure B2. MAIN Program Functional Flow Chart. 
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EEG eG. box 


READ ATRERAET 


PARAMETERS 


READ AERODYNAMIC 
COBFRTCIENTS 





Pra 
DATA? 


ES 







PRINT AERODYNAMIC 
CORFE TCIENTS 


RETURN 


Subroutine DATAIN Function Flow Chart. 


so 





oe U kab 
ASotGn VAbUES TO ,COUNTERS AND FLAGS 


CALL SUBROUTINE FN 











INITIAL 
CONDITIONS 
OUTS DDE 
BOUNDARY? 





"INITIAL ARGUMENTS 
OUTSIDE BOUNDARIES: 














STOPPING 
CONEUTLONS 
om Sh TED. 


F RETURN 


INCREASE STAGE COUNTER 


COMPUTE FIVE ORTHOGONAL 


SECONDARY DIRECTIONS BY 
GR SCuMiDT “PROCESS 





ITERATION DO LOOP 






CALL SUBROUTINE EN TO 
DETERMINE -CRITERIA: FN 
VAGUE, Ad Aa 


Figure B4. Subroutine EXTREM Functional Flow Chart. 
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CALL SUBROUTINE FN TO 
DETERMINE CRITERIA 
FUNCTION AT X + DX 


REDUCE STEP 
S ZE 










OUTSIDE 
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Vet 


Be PARABOLA 


CALI SUBROUTINE: £N. TO 
DETERMINE APPROXIMATE 
MINIMUM (FOPT) 











REDUCE STEP 
SoZ 






OUlo1 DE 
OUNDARY ¢ 










Ee 


SET NEXT ORTHOGONAL SEARCH 
DERECTION 2ND SIEP  sSiZe& 


DETERMINE NEW MAIN DIRECTION 
AND SIEP SiI2Z2b FOR THE NEAL 
STAGE WITH X = XFOPT 


Figure B4. (Continued) 
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Sak 


SEF SEARCH VARIABLES 


CA ER SUBROUTINE CONST 






Ee 


TABLE 
‘Gen 


com 


CALE SUBROUTINE ClLOKUP : 


CALE SUBROUTINE SSEONS 


CALCULATE VALUE OF CRITERIA FUNCTION 
& INCREASE ITERATION COUNTER BY ONE 


. 


RETURN 


Figure B5. Subroutine FN Functional Flow Chart. 
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CALCULATE SPiN 
BON. VARIABLES 
RETURN 


muamce B6. Subroutine CLOKUP Functional Flow Chart. 


START 


DETERMINE INTERPOLATION 
INDICES PROM ALPHA, soBETA 





INTERPOLATE ALPHA & BETA 
PE PENDENT -CORPEPTICIENTS 


INTERPOLATE ALPHA 
DEPENDENT COBPPICliENrTS 


RETURN 


Figure B7. Subroutine CLOKUP Functional Flow Chart. 
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START 


CALCULATE THE 
RED PDUAES FOr THE 
SE UN SEOUATIONS 


RETURN | 


Peaure BS, Subroutine SSEONS Functional Flow Chart. 















wee Ee REPARATION OF THE DATA DECK 

mine LOllowing 1S a guide for preparation of the data 
deck used by computer program SPIN. 

CARDS 1 and 2 — The initial conditions on GAMDOT (rad/ 
Peeve Rift), RDOT(£E/Sec) , PHI (deg), THETA (deg), PSI (deg), 
and DELTAA(deg) on card one and DELTAR(deg), DELTAE(deg), 
THRUST (lbs) and ZDOT(ft/sec) on the second. Data is input 
Ee wolodcind wpOIneraeccondime "towthe format) (7/Fl1iS5/4F11.5). 

CARD tine Peneconeiins bapanerers for Subroutine 
EXTREM. DFMAX, DXMAX, FMAX, LMAX and IW are read in ac- 
Goud nd to thes format (3610, 2,15,12). 

CARD 4 ~— Aircraft identification name (e.g., CONFIGURA- 
THON A) Starting in column one. 

CARDS oeand. 6 — ine macSeand physical characteristics 
of the aircraft, atmospheric density and gravitational 


acceleration values. Card 5 contains: SSAM (mass, Slugs), 
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S (wing area, aay B (wishes pan ,@ttie CBAR (c, ft) ea neces 
(gravitational acceleration, eee Card 6 contains: 

KT (I, slugs 2k YI opt Siege Zo (lo, slugs—ft*), 
EAL ape sities—res) and RHO (atmospheric density, slugs/ 


ie). bata son GCach card are input according to format 


foe 16. 3) . = 

CARDS 7 — 412 — (alternating three card series). These 
cards contain the tabulated aerodynamic coefficients that 
are dependent on both angle of attack and sideslip. The 
first card in the series is as follows: in columns 1&2 
mime COoectficient number (1), columns 384 blank, column 5 
the index on beta (IB), columns 6-10 blank and the remain-—- 
in 70 columns, in 10 column increments, are the coeffi-— 
cients indexed on alpha. The second card of the series is 
divided into eight 10 column increments and the third card 
is divided into four 10 column increments of aerodynamic 
coefficients indexed on alpha. Alpha is indexed 19 times 
for each beta index which corresponds to a range of alpha 
a@OnverOn it Oct: 9A) neGneoet in 5 degree intervals. Beta 1s 
indexed nine times which corresponds to a range of beta 
from —40 to +40 degrees, in 10 degree intervals. 

Lio ii GeacOumOtmtinee Cards, Within a) particular cocik= 
ficient grouping, corresponds to an angle of sideslip equal 
to —40 degrees, where the first coefficient in columns 
11-20 corresponds to an angle of attack equal to 0 degrees, 


columns 21-30 corresponds to an angle of attack equal to 
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5 degrees and so on in 5 degree increments of alpha until 
90 degrees is reached in columns 41-50 of card three. 

The second set of three cards contains coefficients over 
the alpha range for beta equal to ~30 degrees. The remain- 
ing series have coefficients for beta equal to -—20°, ~—10°, 
Umer lO°, +20°, +30° ane aoe The next grouping of 27 
cards (nine three card series) contains another aero- 
dynamic coefficient over the range of alpha and beta and 
Semen in 2/ card groupings until all 15 alpha and beta 
dependent coefficients have been prepared. Format for the 
mmeee card Series 1S (12,2X,5£1,5x%,7F10.7/8F10.7/74F10.7). 

The input Sequencing of groupings 1S very important as 
the machine computation is dependent upon the assignment 
of a particular tier in the three~—dimensional COEFAB array 
mowENe proper aerodynamic coefficient. Table Bl gives the 
proper assignment. 

CARDS 413-429 — (alternating three card series). These 
cards contain the tabulated aerodynamic coefficients that 
are dependent a alpha only. The first card in the series 
is as follows: columns 1&2 the coefficient number (I), 
columns 3-10 blank and the remaining 70 columns, in 10 coli- 
umn increments, are the coefficients indexed on alpha. The 
second and third cards of the series are the same as de~- 
scribed above for cards 6-411. Again alpha is indexed 19 
Elimes "COmresponding EO a range ©f alpha from 0 £o 90° in 


5° increments. Columns 11-20 of card one corresponds to 


5S 








the coefficient for alpha equal to zero degrees, columns 
Zi—30 for alpha equal to 5°, and so on in 5° increments 
until alpha equal to 90° is reached in columns 41-50 of 
meard three. Format for these three card series is (12,8x, 
PeetOe//SFl0.7/4F10.7). As before, the input sequence of 
mre coerficient groupings is very important for proper 
interpretation by the program. The assignments are given 


in Table Bl. 
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Coefficient Assignment 


Storage Coectruclent 
Be@erricient bhoecation Preinw sOut 
7 


1 
COnnAs (el 7) 
COEPABR (72, ,_) 
COBRAB ( 3, _,_) 
COEFAB ( 4, ,_) 
C@nuAss ( 5, 5) 
COHEAB (267. 47) 
COEFAB ( 7, ,_) 
CORE ABI Oy. 71.) 
COEEABaK Io +)? 
COLT AB MLO.) 
COEFAB (11,7; _) 
COBH ABERGL2Z 7) ) 
COEGEAP S13 c78 9) 
COEFAB (14, _,_) 
COU EA Bea Gio a.) 
COErZ@ (ltr) 
CCnt am (62) 
COEFA ( 3, _) 
COFFA ( 4, ) 
COEEAS (oop) 
COEFA ( 6, ) 
COERA MS (/7ae) 
COEE AS A. 87) 
CORDA ter (9738) 
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Location 


CR( 7) 
CR( 8) 
CR( 9) 
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CR(14) 
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CR(16) 
CR(17) 
CR(18) 
CR(19) 
CR (20) 
CR(21) 
CR (22) 
CR (23) 


CR(24) 
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G. PROGRAMMING CONSIDERATIONS 

Program SPIN requires 66K bytes of storage. When buffer 
Space 1S taken into account, a request for 80K bytes should 
be made when submitting the job for processing. 

The time requirement 1S variable and is very dependent 
on the closeness of the eel conditions to the steady 
Spin condition. Convergence has been obtained with a total 
CPU time of a little over one minute to as high as 15 min- 
utes; this time includes about 25 seconds for the compile 
and link steps. 

The optimization scheme loop time is very rapid. The 
time required to complete one stage is 0.2 second, which 
meemires 18 criteria function evaluations. The criteria 
mumetlon evaluation occurs at the rate of one per .0111 
second and this accomplishes, as Figure Bl indicates, the 
evaluation of the equation constants, interpolation of the 
24 aerodynamic coefficients, and the evaluation of the 
equations of motion. By looking at the listing for each 
of the aurea nee involved, one can better appreciate the 
speed with which the program operates. 

The variable time requirement of the program can create 
pmpotoplLem to the researcher. More often than not, the 
turn around time for any program is determined by the time 
and storage requested. As a general guideline to the use 
of program SPIN, the following is suggested: 

a. Set LMAX (the maximum number of stages) to a value 


Stat, O50. 
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Request 80K storage and a four minute oe time. 
Chee eresiltsS Omeeine run, If the value of the 
Sclecmanruncelton 1S qreater than Ona guess 
another set of initial conditions and rerun. Re- 
Snesominaveate that 1£ the value of the criteria 
minction hasn't dropped below the above value in 
the four minutes alloted that an increase of time 
will not yield a steady spin condition. The pron 
Gaaterma this case 1s convergingeon a local mini- 
mum which is not indicative of a steady spin. If 
eniesvalue Of the criteria function is less than 
107, rerun the program with LMAX set at 2,500 
and a requested time of nine minutes if greater 


accuracy is desired. 
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eee Nie 
AIRCRAFT DATA 


Hieerollowing aircraft data was used in the study of 
steady spin conditions. Table Cl tabulates the mass and 
dimensional characteristics of the two aircraft and the 


following pages list the aerodynamic coefficients for each 


Semraguration. 
fap cil. Mass and Dimensional Characteristics 
Characteristics Gomi ilouraclonss Conitlvquration - 
m (slugs) 1554500 Dy eee 
S ee) D220 0 695.00 
ip (LE) 63.00 38.00 
er (Et) 9.04 232796 
I, (slugs—ft*) 53100.00 13600.00 
I, (slugs-£t*) ~ | 299000.00 128000.00 
I (slugs-ft*) S30750.00 138000.00 
I, (Slugs-£t*) 12480.00 4340.00 
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